﻿class verlet
{
private:
	caja *box;
	//Paso de tiempo Dt
	double Dt;
	//Dt2
	double Dt2;
	//Dt/2
	double Dth;
	//Dt2/0.5
	double Dt2h;
	//Uso del termostato
	bool termo;

public:

	verlet(caja *bx, double xDt, bool xtermo)
	{
		box=bx;
		Dt=xDt;
		Dt2=Dt*Dt;
		Dth=Dt*0.5;
		Dt2h=0.5*Dt*Dt;
		termo=xtermo;
	}

	void  LeapFrogLinkCellList(void);
	void  verletLinkCellList(void);

};

void  verlet::LeapFrogLinkCellList(void)
{
	double dv[3];
	double parke=0.;
	double posR;
	int mc[3], c;

	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;



	for(int i=0;i<box->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			box->part[i].pos[j]+=(box->part[i].vel[j])*Dt+(box->part[i].force[j])*Dt2h;
			//Calculamos la velocidad en t+Dt/2
			if(termo)
			{
				//Empleamos un reescalado de las velocidades como termostato
				box->part[i].vel[j]=(box->part[i].vel[j])*sqrt(box->temp0/box->temp)+Dth*box->part[i].force[j];
			}
			else
			{
				//No utilizamos termostato
				box->part[i].vel[j]+=Dth*box->part[i].force[j];
			}
		}
		box->boundary(i);
	}
	//Calculamos las fuerzas, epotencial y el virial
	interaccionLinkCellList(box);

	//Reset the head of the link-cell list
	for (int i=0;i<box->ncell3;i++)
	{
		box->head[i] = -1;
	}
	//Velocidad en t+Dt
	for(int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]+=(box->part[i].force[j])*Dth;
			//Actualizamos la componenete de la energia cinetica
			parke+=SQR(box->part[i].vel[j]);
			//Calculamos la deriva media de las velocidades de la caja
			dv[j]+=box->part[i].vel[j];
			//Coordenadas de la nueva celda en la que se encuentra la particula
			mc[j] = int((box->part[i].pos[j]+0.5*box->lbox[0])/box->lgrid);
			box->part[i].celdaxyz[j]=mc[j];
		}
		///Indice de la celda donde se encuentra la particula
		c = int(mc[0]*box->ncell2+mc[1]*box->ncell+mc[2]);
		box->part[i].celda=c;
		///Link to the previous occupant (or EMPTY if you're the 1st)
		box->link[i] = box->head[c];
		///The last one goes to the header"""
		box->head[c] = i;
	}

	///Calculamos la velocidad media de la caja
	for(int i=0;i<3;i++)
		{box->vel[i]=dv[i]/double(box->npart);}
	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/double(box->npart);
	//Calculamos la temperatura de la caja
	box->temp=parke/3.0/double(box->npart);
	///Calculamos la presion de la caja
	box->presion=box->rho*box->temp+box->virial/box->vol; //IMPORTANTE-> Falta introducir la correccion del cutoff a la presion
}

void  verlet::verletLinkCellList(void)
{
	//double posNew[3];
	double dr;
	double dv[3];
	double parke=0.;
	double posR;
	int nper;
	double anper;
	int mc[3], c;

	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;

	//Reset the head of the link-cell list
	for (int i=0;i<box->ncell3;i++)
	{
		box->head[i] = -1;
	}

	for(int i=0;i<box->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			posR=(box->part[i].pos[j])+(box->part[i].vel[j])*Dt+(box->part[i].force[j])*Dt*Dt/2.0;
			//Condiciones de contorno periodicas
			nper = int(dabs(posR)/(0.5*box->lbox[j])); //si esta dentro de la caja 0; si esta fuera 1
			anper = sign((double(nper)*box->lbox[j]),posR); //0 dentro de caja, 1 � -1 si esta fuera
			posR-=anper;

			//posNew[j]=(bx->part[i].pos[j])+(bx->part[i].vel[j])*Dt+(bx->part[i].force[j])*Dt*Dt/2.0;
			//dr=posNew[j]-(bx->part[i].posOld[j]);
			dr=posR-(box->part[i].posOld[j]);
			dr-=int(dr/(0.5*(box->lbox[j])))*(box->lbox[j]);
			box->part[i].vel[j]=((dr)/(2.0*Dt));
			box->part[i].posOld[j]=box->part[i].pos[j];
			//bx->part[i].pos[j]=posNew[j];
			box->part[i].pos[j]=posR;
		}
		box->boundary(i);
	}

	///Uso del termostato?
	if(termo)
	{
		vscale2(box);
	}
	for(int i=0;i<box->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			parke+=SQR(box->part[i].vel[j]);
			dv[j]+=box->part[i].vel[j];
		}
	

		///Creamos la link-cell list
		for(int j=0;j<3;j++)
		{
			mc[j] = int((box->part[i].pos[j]+0.5*box->lbox[0])/box->lgrid);
		}
		c = int(mc[0]*box->ncell2+mc[1]*box->ncell+mc[2]);
		///Link to the previous occupant (or EMPTY if you're the 1st)
		box->link[i] = box->head[c];
		///The last one goes to the header"""
		box->head[c] = i;
		
	}
	for(int i=0;i<3;i++)
		{box->vel[i]=dv[i]/double(box->npart);}
	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/box->npart;
}

